########################################################################
# PLS 206 — Week 4 live coding, part 1: curves
# Grey Monroe, Fall 2026
########################################################################

# install.packages(c("ggplot2", "mgcv"))
library(ggplot2)
library(splines)
library(mgcv)

# fitted values with a 95% CI, for any lm() or gam()
fit_ci <- function(model, newdata) {
  p <- predict(model, newdata = newdata, se.fit = TRUE)
  data.frame(newdata, fit = as.numeric(p$fit),
             lwr = as.numeric(p$fit - 1.96 * p$se.fit),
             upr = as.numeric(p$fit + 1.96 * p$se.fit))
}

########################################################################
# 1. The data
########################################################################

# Wheat seedlings, one pot each, 14 days at 8-42 degrees C.

seedlings <- read.csv("https://greymonroe.github.io/PLS_206/data/seedling_temp.csv")
std <- subset(seedlings, variety == "standard")
head(std)
nrow(std)

ggplot(std, aes(temp, mass_mg)) +
  geom_point() +
  theme_classic()

########################################################################
# 2. A straight line and its residuals
########################################################################

m_lin <- lm(mass_mg ~ temp, data = std)
summary(m_lin)      # slope p = 0.59: the line sees nothing

grid <- data.frame(temp = seq(8.3, 42, by = 0.05))
ci_lin <- fit_ci(m_lin, grid)
ggplot(std, aes(temp, mass_mg)) +
  geom_point() +
  geom_ribbon(data = ci_lin, aes(temp, ymin = lwr, ymax = upr), inherit.aes = FALSE, alpha = 0.3) +
  geom_line(data = ci_lin, aes(temp, fit), colour = "blue") +
  theme_classic()

# residuals: always look
plot(std$temp, resid(m_lin))
abline(h = 0, lty = 2)

########################################################################
# 3. A quadratic: one new column
########################################################################

m_quad <- lm(mass_mg ~ temp + I(temp^2), data = std)
summary(m_quad)     # temp^2 negative: bends down

ci_quad <- fit_ci(m_quad, grid)
ggplot(std, aes(temp, mass_mg)) +
  geom_point() +
  geom_ribbon(data = ci_quad, aes(temp, ymin = lwr, ymax = upr), inherit.aes = FALSE, alpha = 0.3) +
  geom_line(data = ci_quad, aes(temp, fit), colour = "blue") +
  theme_classic()

# the vertex: estimated optimum
b <- coef(m_quad)
-b[2] / (2 * b[3])  # 24.6 degrees; the true optimum is 26 (we simulated it to be so)

plot(std$temp, resid(m_quad))
abline(h = 0, lty = 2)

########################################################################
# 4. Higher polynomials
########################################################################

m_p10 <- lm(mass_mg ~ poly(temp, 10, raw = TRUE), data = std)
summary(m_p10)$r.squared

# beyond the data
wide <- data.frame(temp = seq(4, 46, by = 0.1))
ci_p10 <- fit_ci(m_p10, wide)
ggplot(std, aes(temp, mass_mg)) +
  geom_point() +
  geom_ribbon(data = ci_p10, aes(temp, ymin = lwr, ymax = upr), inherit.aes = FALSE, alpha = 0.3) +
  geom_line(data = ci_p10, aes(temp, fit), colour = "red") +
  coord_cartesian(ylim = c(-40, 120)) +
  theme_classic()

########################################################################
# 5. A spline: bumps as columns
########################################################################

X <- bs(std$temp, df = 6)
round(head(X), 2)   # one row per pot, one column per bump

m_bs <- lm(mass_mg ~ bs(temp, df = 6), data = std)
round(coef(m_bs), 1)   # one weight per bump

ci_bs <- fit_ci(m_bs, grid)
ggplot(std, aes(temp, mass_mg)) +
  geom_point() +
  geom_ribbon(data = ci_bs, aes(temp, ymin = lwr, ymax = upr), inherit.aes = FALSE, alpha = 0.3) +
  geom_line(data = ci_bs, aes(temp, fit), colour = "blue") +
  theme_classic()

# a natural spline: straight lines beyond the data, and it drops into any model formula
m_ns <- lm(mass_mg ~ ns(temp, df = 6), data = std)
wide <- data.frame(temp = seq(0, 50, length.out = 300))
ci_ns <- fit_ci(m_ns, wide)
ggplot(std, aes(temp, mass_mg)) +
  geom_point() +
  geom_ribbon(data = ci_ns, aes(temp, ymin = lwr, ymax = upr), inherit.aes = FALSE, alpha = 0.3) +
  geom_line(data = ci_ns, aes(temp, fit), colour = "blue") +
  theme_classic()

########################################################################
# 6. Cycles: sin and cos columns
########################################################################

soil <- read.csv("https://greymonroe.github.io/PLS_206/data/soil_season.csv")

f1 <- lm(flux ~ sin(2 * pi * day / 365) + cos(2 * pi * day / 365), data = soil)
f2 <- lm(flux ~ sin(2 * pi * day / 365) + cos(2 * pi * day / 365) +
                sin(4 * pi * day / 365) + cos(4 * pi * day / 365), data = soil)
summary(f1)$r.squared   # 0.60
summary(f2)$r.squared   # 0.73

days <- data.frame(day = 1:365)
ci_f2 <- fit_ci(f2, days)
ggplot(soil, aes(day, flux)) +
  geom_point(colour = "grey60") +
  geom_ribbon(data = ci_f2, aes(day, ymin = lwr, ymax = upr), inherit.aes = FALSE, alpha = 0.3) +
  geom_line(data = ci_f2, aes(day, fit), colour = "blue") +
  theme_classic()

########################################################################
# 7. A first GAM
########################################################################

m_gam <- gam(mass_mg ~ s(temp), data = std, method = "REML")
summary(m_gam)      # edf 6.87

ci_gam <- fit_ci(m_gam, grid)
ci_gam$temp[which.max(ci_gam$fit)]   # fitted optimum: 26.25 degrees

ggplot(std, aes(temp, mass_mg)) +
  geom_point() +
  geom_ribbon(data = ci_gam, aes(temp, ymin = lwr, ymax = upr), inherit.aes = FALSE, alpha = 0.3) +
  geom_line(data = ci_gam, aes(temp, fit), colour = "blue") +
  theme_classic()

# k is a ceiling, not the wiggliness
summary(gam(mass_mg ~ s(temp, k = 5), data = std, method = "REML"))$edf    # 3.95: k too small
summary(gam(mass_mg ~ s(temp, k = 20), data = std, method = "REML"))$edf   # 7.75